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Abstract 

We consider sensor array imaging for simultaneous noise blended sources. We study a mi- 
gration imaging functional and we analyze its sensitivity to singular perturbations of the speed 
of propagation of the medium. We consider two kinds of random sources: randomly delayed 
pulses and stationary random processes, and three possible kinds of perturbations. Using high 
frequency analysis we prove the statistical stability (with respect to the realization of the noise 
blending) of the scheme and obtain quantitative results on the image contrast provided by the 
imaging functional, which strongly depends on the type of perturbations. 
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1 Introduction 

An amazing fact in the analysis of imaging functionals, that has been recently pointed out and 
is currently under investigation, is that specific kinds of noise can improve the image quality or 
drastically reduce computational costs associated to the evaluation of the imaging functional. A 
strong motivation is the observation that time reversal refocusing is enhanced when the medium 
is randomly scattering. A time reversal experiment is based on the use of a special device called 
a time reversal mirror (TRM), which is an array of transducers that can be used as receivers and 
transmitters. A typical time reversal experiment consists in two steps. In a first step, a point source 
emits a short pulse that propagates through the medium and is recorded by the TRM used as an 
array of receivers. In a second step, the recorded signals are time reversed and reemitted by the TRM 
used as an array of transmitters. The waves then refocus on the original point source location. The 
striking observation is that refocusing is enhanced when the medium is randomly heterogeneous and 
scattering compared to the situation in which the medium is homogeneous. Moreover, the refocused 
pulse is statistically stable in the sense that it does not depend on the realization of the random 
medium, but only on the statistical distributions of the fluctuations of the random medium [FGPS]. 

In the context of sensor array imaging, a similar technique is employed. The typical experiment 
still consists in two steps. The first step is the experimental data acquisition: a point source emits 
a wave into the medium, the wave is reflected by the singularities in the propagation speed of the 
medium and is recorded by an array of receivers. The second step is a numerical processing of 
the recorded data: the recorded signals are time reversed and resent in a numerical simulation into 
a model medium to locate the singularities in the propagation speed. However, in contrast with 
physical time reversal, the fictitious medium employed in the imaging process cannot in general 
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capture complex heterogeneities of the original medium. Therefore, research has focused on other 
approaches to improve the imaging process, such as the use of random sources [DFGS12], [GP10], 
[HSH08], [SWH11], [VB11], [WNT12]. 

The classical approach to the imaging problem consists in performing a large number of exper- 
iments sounding each time a different source. For each experiment, the signal recorded at each 
receiver is stored and time reversed. This produced the data matrix. For each experiment, the 
time reversed data are reemitted (numerically) into the model medium in a new simulation, and the 
images obtained by each simulation are stacked. While this method provides a very good "image" of 
the medium, it involves gathering, storing and processing huge amounts of data [Bc09], [MDB11]. 

In the noise blending approach noisy sources are used and they are all sounded simultaneously 
in one experiment. In this case, the time reversed recorded signals from the physical experiment are 
stored in a data vector and resent simultaneously into the model medium in a single simulation. This 
approach allows for considerable savings in both the data gathering, storing and processing stages. 
However, special care must be put in the choice of the noisy sources in order to ensure that the cross 
talk terms are very small and do not compromise image quality. This technique can be successfully 
applied also to the time reversal approach and its analysis bears similarities with techniques for 
passive imaging, which exploits ambient noise sources to recover travel times from correlations in 
between recordings at different stations [GP09]. Applications of these different techniques are being 
investigated in different fields, from seismology [SCS06] [LMD06] [SCSR05] [GSB08], to volcano 
monitoring [SRG06] [BSC08] [BSC07], to petroleum prospecting [CGH06] and medicine [FCDOO]. 

As remarked, the crucial step for the noise blending approach lies in the choice of the noisy 
sources. In [DFGS12] it was suggested the use of stationary random Gaussian sources or of randomly 
delayed source pulses: for these choices of random sources, fourth moment computations show that 
the algorithm is statistically stable. In the present work we develop some further investigations on 
this setting and show how to obtain quantitative results on image quality and statistical stability 
of this algorithm in the high frequency regime, when the goal is to image singular perturbations in 
the speed of propagation. The result strongly depends on the type of perturbation. We therefore 
consider three types of perturbations, supported respectively on small balls, thin tubes and thin 
discs. With a slight abuse of notation we will refer to them as point, line and surface singularities, 
as they can be thought of as approximations of singular perturbations of the velocity of propagation 
supported on subspaces of lower dimension. 

For each kind of perturbation we first analyze the average image contrast seen between the 
center of the perturbation and a point far from it: this will provide an hint on the level of difficulty 
to correctly image these perturbations. An even more interesting result follows: it concerns the 
quantitative analysis of the statistical stability of this functional, providing the typical contrast seen 
for the three perturbations. The question of stability of the imaging functional has already been 
addressed in [DFGS12]: no quantitative analysis was carried out there, but it was shown that a 
condition for the statistical stability is that the recording time interval T must be large. With a 
careful analysis of fluctuations produced by the random sources we show that, in the high frequency 
regime, the typical contrast is actually much better than just of order y/T; but again the exact order 
of amplitude of the contrast depends on the shape of the perturbation. Point singularities are easy 
to observe, while surface type perturbations are the hardest to locate. 

The paper is organized as follows. A short presentation of the model and the imaging functional 
used will complete this introductory section. In section 2 we analyze the average (with respect to 
the realization of the random time delays used in the blending process) sensitivity of the method to 
the three types of singularities. In section 3 we study fluctuations due to the stochastic nature of the 
result, and from the analysis of the typical behavior we obtain conditions ensuring the possibility to 
accurately image the perturbations. Finally, all the results obtained are collected and discussed in 
section 4. 
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Notation. We use boldfaced characters to denote vectors: for example but x = (x, y, z) £ R 3 . 

B r will denote the ball of radius r in R d for d equals 2 or 3: B r = {x £ M. d \ \x\ < r}. The notation 
A = 0(e) means that the quantity A is exactly of order e; in order to say that it is of order e or 
smaller we will write A < 0(e). 

1.1 The wave equation 

We consider the solution u of the wave equation in a three-dimensional inhomogcneous medium 

1 d 2 u , 

^)W- AxU = nit ' x) ' (L1) 

where c(x) is the velocity of propagation of waves in the medium and x = (x, y, z) £ R 3 . We rewrite 
the velocity in the form 

c~ 2 (x) = Cq 2 (x) + 5c- 2 (x) , 

where co(x) is the known smooth background velocity (for simplicity we assume it to be constant) 
and 6c~ 2 (x) is the velocity perturbation that we want to estimate, whose spatial support is contained 
in some domain C i 3 . We will detail these perturbations in the following. To simplify the 
geometry of the model, we shall take ft = Br to be a ball with a large radius R. We also assume 
to have N B point sources located at points (y s )s=i,...,N 3 laying on the surface dB R . They can either 
emit (almost) simultaneously the same short pulse waveform, but randomly delayed in time, or 
independent stationary random signals. Following [DFGS12], we will refer to the first case as (noise) 
blended sources and to the second as stationary random sources. 
For the first case, the source term n(t,x) is of the form 

N B 

n(t,x) = f(t - t s )S{x - y s ) . 

8=1 

The (short) pulse function (/(t))tgR is deterministic. Its carrier frequency is u>o and its bandwidth 
is b. The time delays (T s ) s= i i ... i jv s are zero-mean independent and identically distributed random 
variables with probability density function p T (t). 
For the second case, the source term is given by 

N s 

n(t,x) = ^2n s (t)S(x - y s ) , 

8=1 

where the random functions (n s (i)) tg R, s = 1, . . . , N s , are independent, zero-mean, stationary Gaus- 
sian processes with autocorrelation function 

(n s (t x )n s -(t 2 )) =8 ss ,F(t 2 -t 1 ). 



The direct and inverse problems can be formulated in terms of the background Green's function, 
the fundamental solution of the wave equation (1.1). For an homogeneous medium with constant 
background velocity cq, in the Fourier domain the Green's function is given by 

G(u,x 1 ,x 2 ) = — ; — r exp (i— \x 1 - x 2 \) ■ 

4ir\xi — x 2 \ V co / 

Here, the Fourier transform of a function f(t) is defined by 

/>) = J f(t)e iut dt. 
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1.2 Direct and inverse problems 

We introduce here the scattering operator, that is, the mapping from velocity perturbations to the 
data, in the Born approximation [BCS]. 

We assume that sources (located at points y s , s = 1, ...,N S ) are disposed on the surface of the 
ball Br containing the perturbations, and are dense enough (ideally, closer than half of the central 
wavelength) so that a continuum approximation can be used. Signals are observed at the passive 
sensor array (x r ) r =i,...,N r for some large time interval [— T/2, T/2]. For noise blended sources, the 
recording time T should be much larger than the typical travel time, so as to guarantee that the 
backscattered signals are completely recorded. For stationary random sources it must be taken much 
larger that the inverse of the bandwidth of the noise sources (i.e. much larger than the decoherence 
time). 

The recorded data consists of the vector d(t) = (d(t, x r )) r —i ... jv r of the signals recorded by x r 
for t € [—T/2, T/2]. These data arc modeled by the scattering operator T : (8c~ 2 (x)} xeB —> 

( d (*)) ie[ -T/2,T/2]' Wh(3re 



(FSc 2 )(t,x r ) = / Q(t,x r ,x)Sc 2 (x)dx. 
Jn 



Q(t, x r , x) = -^ JffG(tux r , x)G(t 2 ,x, y)n(t -t x - t 2 , y) dh dt 2 dy . (1.2) 
In the Fourier domain 

Q(u), x r , x) =u? J G(u, x r , x)G(u, x, y)h(u), y) dy . (1.3) 

This is the formulation of the direct problem: the expression of the data set in terms of the velocity 
perturbation. 

The imaging problem (inverse problem) aims at inverting the map J- in order to reconstruct the 
velocity perturbation (<5c~ 2 (a;)) xEBr from the data set (d(i)) tg j_ T y 2 T , 2 y The usual (least-square) 
approach would consist in applying the operator (J 7 * \F)~ l T* to the data set d, where the adjoint of 
the scattering operator T is given by 

N r I. 

{T*d){x) = ( 2 Q{t,x r ,x)d{t,x r )dt. (1.4) 

r=l _ T 

However, the full least-square inversion is in practice too complicated and the normal operator 
T*J- is usually dropped in the inversion process. In [DFGS12] it was shown that for T large the 
normal operator is statistically stable (i.e. its fluctuations are smaller than its expectation) and 
that its statistical average is close to the identity operator (more precisely, the kernel (T* T)(x,x') 
concentrates near the diagonal x = x'), proving that this procedure can indeed provide a reasonable 
estimate of the velocity perturbation. 

In the following sections we perform a detailed analysis of this imaging functional in the high 
frequency regime, obtaining quantitative results on the statistical stability of the method for different 
types of perturbations. 

Since the main advantage of our approach lays in the drastic reduction in storage and computa- 
tional costs, let us stress that the computation of the adjoint operator J-* can be done quite easily. 
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Remark 1. The adjoint operator T* can be computed solving only two wave equations as follows. 
First, compute the wave u(t, x) emitted by the original source, which is to say solve the wave equation 
with source n(t, x) and background velocity cq: 

u(t,x) = [j G(t 1 ,x,y)n(t-t 1 ,y)dydti . 

Second, compute the anti causal solution v(t,x) of the wave equation with source term 
J2 r S{y ~ x r )dfd(t,x r ): 

N r . N r . 

v(t, x) = / G(t 2 - t,x,x r ) dt 2 d(t 2 ,x r )dt 2 = ^ / G(t 2 , x, x r ) dfd(t 2 + t, x r ) dt 2 . 

r=l*' r=l*' 

Correlating the two wave solutions at a point x in the search window produces the imaging functional 

z W = -/„„,*)„(,, *)d, 

= X] J Jjj G(h,x,x r )G(t 1 ,x,y)n(t-ti-t2,y)d(t,x r )dt 1 dt 2 dy dt 
and using the definition (1.2) we obtain 

l{x) = {F*d){x) . 
1.3 Analysis of the imaging functional 

Let us start by the analysis of a kernel which will appear in the imaging functional we have to study. 

In our model, sources and receivers are located on the surface of the ball Br containing the 
perturbations and are dense enough so that a continuum approximation can be used. We can 
therefore rewrite the kernel 

fC u (x,y) := y2&(u,x r ,x)G(u;,x r ,y) ~ / p(x r )G(u>, x r , x)G(u, x r , y) da(x r ) , 

where p{x r ) is the surface density of receivers or sources. We also assume that perturbations are 
located near the center of the ball; this means that the medium is homogeneous outside of a ball B r 
of radius r <gi R. Then, we have the approximate identity [GP09, Proposition 4.3] 

— / G(oj, x r , x)G(w, x r , y) da(x r ) ~ 6{u, x, y) - G(u, x, y) , 
c o JdB R 

which follows from Green's identity and the Sommerfcld radiation condition. This result can be 
viewed as a version of the Helmholtz-Kirchhoff integral theorem. Using the function sinc(x) = 
sin(x)/x, the right hand side of the above equation can be rewritten as 

2i 3( G(uj, x, y) | = — - — sincf — \x — y 
V / 4?r c Vc 

and assuming that receivers and sources have constant density on the surface dBn we obtain 

£u,(x,y) = —sincf— \x - y\) . 

47T V Cq / 
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We return to the analysis of the imaging functional, which estimates the velocity perturbations, 
given by (1.4). Using equation (1.3) and the definition of the vector d(t) of recorded signals, we can 
write the imaging functional as 



Sc 2 {x) = (F*<£)(x) = — f I uj 2 G(u, x, x r )G(tj, x, y)h{uj, y)d(uj, x r ) dydui 

= [ F*F(x,x r )Sc- 2 (x')dx' = [ JC(x,x r )6c- 2 (x')dx' , (1.5) 

J Br J Br 

where 8cr 2 (x) is the estimated velocity perturbation. We will use (•) to denote the statistical average 
with respect to the distribution of the random time delays or the stationary random source signals. 
The statistical average of the kernel K, (which is not the same kernel as the K u used earlier) is given 
by 

(K(x,x')) = (T*T){x,x') 

= — uj 4 \f(uj)\ 2 ]T \6(cj,x,y s )G(uj,x',y s )] ]T [G(w,x r ,x)G(u,x ri x')\ du 

J s=l r=l 

= ^j" A \m\ 2 ^ c 2 (^\x-x>\)d. 

for noise blended sources, and by 

(JC(x, x')) = -^r ( w 4 F(w)sinc 2 ( — \x - x'\] dw 
2°it 6 J Vco / 

for the stationary random sources, see [DFGS12]. 

We stress that the focus of this paper is on the high frequency analysis of this functional. We 
will analyze its performances in localizing singularities in the background velocity of propagation, in 
the high frequency regime 77 <gc 1 , where 

wo _ 1 
Co V ' 

Note that the wavelength is 2-kt). 



2 Expected contrast of the estimated perturbation 

Even if the estimated perturbation Sc~ 2 (x) provided by the imaging functional does not have the 
exact shape of the real perturbation, due to the different approximations used, it still shows a peak 
on the actual location of the velocity perturbation. We study in this section the expected (average 
in the statistical sense) contrast seen in the estimated velocity perturbation between the location of 
the real perturbation and points far from it. 

The expected estimated perturbation for noise blended sources is given by 

(5c- 2 {x)) = -L y Jifow^l/Mfdw , 

with 

1 1 (x,uj)= / sinc 2 (— \x - x'\]Sc- 2 (x')dx' . 
Jb r Vc o > 
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When the bandwidth of the source pulse is smaller than the central frequency, it is enough to study 
the behavior of the spatial integral X\ (x) at the central frequency. For simplicity we drop the 
dependence on uj. It turns out that T ■ T\(x) is the quantity one has to study also in the case of 
stationary random sources. All the computations presented in the remaining part of this section 
are performed in the noise blended case; to obtain the different orders of amplitude provided by the 
imaging functional for stationary random sources it suffices to multiply the results by T. 

To simplify the presentation of some computations, we assume that the support of the pertur- 
bation is centered on the center of the ball Br on the surface of which are located the sources and 
receivers, and choose this center as the origin of our system of coordinates. 

For the purposes of this section it would not be necessary to distinguish between the scale of the 
wavelength of the signals (ry) and that of the size of the perturbation [e). However, this distinction 
turns out to be crucial in the analysis of fluctuations carried out in section 3. 

2.1 Point singularities 

Let us start by considering a point singularity, that is to say, a singularity whose support has a very 
small diameter. We model it with a perturbation of the velocity of propagation that is supported on 
a ball of radius e: 

6c- 2 (x) = al {Be} (x). (2.1) 

To simplify computations, we take a to be constant. Then, one observes that it only enters formulas 
as a multiplicative constant (squared in section 3): since it is of no relevance to our scopes, we set it 
equal to 1, also for the other types of perturbations. 
Changing variables, we have: 

X 1 (x)= I sine 2 (Ix-x'l/ 7]) 5c- 2 (x' ) dx' = I sine 2 (\x - x'\/i]) dx' 

= e 3 / sinc 2 ( -\x/e - x'\ ) dx' . 
Jb 1 v ?7 1 

At the center of the perturbation, x = 0, we have 

Xi(0) = s 3 J sinc 2 (-|:E'|) dx' = Aire 3 J ^ sin 2 (-p) dp = 2i:erf(l - S mc(2e/r,)). 

For e ~ 7] we have that Xi(0) = 0(erj 2 ). But if e <§; 77, the order becomes 0(e 3 ). 

For I £U I = O(l), what gives the order of amplitude of X\(x) is the decay of the sinc(x) function, 
which goes approximately as l/x. We have that 

l^xj-e 3 [ r] 2 \x -ex'\- 2 sine 2 (\x\/rj) dx' <— e 3 ry 2 . 

Remark that the bound found is sharp, since there are x for which Ii is exactly of order e 3 rj 2 . 

We see that the difference in amplitude between the centre of the perturbation and a point far 
from it is significant: it is of order r\~ 2 when e <ti r], e~ 2 if e ~ r\. 

2.2 Line singularities 

We consider now line-type singularities, which is to say an almost one-dimensional perturbation of 
the velocity of propagation. We model it by a perturbation supported on a cylinder of radius s: 

5cr 2 {x) = 1 {Ce} (x), (7 £ = B £ x[-l,l]cl 2 xl. (2.2) 
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We have 



li(x) = { smc 2 (\x-x'\/rj)5c- 2 (x')dx' 
Jb„ 

{-y/{x - ex') 2 + (y - ey') 2 + (z - z') 2 ) da;' . 



= s 2 / / sinc 2 ( 

J -l J By 

At x = this term reduces to 



e 2 J dz 1 J sine 2 Q-^e 2 {x' 2 + y' 2 ) + z' 2 ) Ax' Ay' , 
which is rapidly oscillating in z' . In the integral in z' we can therefore approximate sin 2 by its mean: 
li(0) ~ 2 7 re 2 ^ 1 Q J ?/ 2 (e 2 r 2 + z' 2 ) _1 dz^rdr = 27re 2 / ^ arctan f i-) r dr 



2-KEif I arctan I — \ Ar = 0(e rf 



er / 



Far from the perturbation the integral is of order e 2 r/ 2 . For example, for x 2 + y 2 = C — 0(1) we 
have 

I. (») £ j£ 2 /' ^ - ^') 2 + (if - + - *')') "' <W 

Therefore, the difference in amplitude seen between the centre of the line perturbation and a 
point far from it is of order e _1 . 

2.3 Plane singularities 

Let us consider now singularities that are approximately two-dimensional: we call them surface- 
type singularities and model them by a perturbation of the velocity which is supported on a disc of 
thickness e: 

Sc~ 2 (x) = t {De} (x), D s = [-£,£] x B 1 c M x R 2 . (2.3) 

With the notation introduced above we have 

Xi(x) = £ J J sine 2 y/Jx - ex 1 ) 2 + (y - y') 2 + (z - z') 2 ) dx' . 



At x = 



Xi(0) = 27re y dx'^ sinc 2 ^VfV+^) rdr 



?re / dx' I n 2 (e 2 x' 2 + r 2 ) r Ar = neif [ ^ In ( 1+ -^-r ) d.x' 



l Jo x J-i 2 



= 0( £?? 2 |ln( £ )|) 
because the sin 2 in the first line is rapidly oscillating in r. 
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For x = 0(1) we have instead that 

Zi{x)<s J J V 2 ((x - ex') 2 + (y - y') 2 + (z - z'f ) 

r-l 



2-rterf 



dr = 0{erf 



, x + r z 

For this type of singularities, the difference in amplitude seen between the centre of the perturbation 
and a point far from it is very weak: it is only of order |Zn(£)|. This already hints to the fact that, 
with the imaging functional we consider, surface-type singularities are harder to locate than the 
other types. 



3 Fluctuations 

We have obtained in the previous section the average contrast of the imaged perturbation. We want 
now to find confidence intervals for the typical contrast observed during the experiment. To do so, 
we must analyze the fluctuations in the result. They are given by the standard deviation of the 
estimated perturbation <5c -2 , which at a point x is given by 



S(x) = ( \Sc- 2 (x) - (6c- 2 (x))\ 



(3.1) 

We need to compute the standard deviation at the location of the perturbation and at points far 
from it, and compare them with the expected amplitude of the estimated perturbation. Using (1.5) 
to write explicitly (3.1) we get 

S(x) = ' 



(K,(x,x') - (K.(x,x')))5c- 2 (x')dx' 
i(fC(x,x') - (fC(x,x')^ (K{x,x") - {K(x,x")\ \ 5c- 2 {x')5c- 2 {x")dx' dx" 
Cov(/C(cc, x'),K{x, x")) 5c- 2 {x')5c- 2 {x") dx' dx" 



In [DFGS12] a formula was obtained for the variance of the kernel K,. It is possible to carry out the 
same computations for the covariance: for noise blended sources one obtains 



27rT T Covl K >C(x,x'),JC(x,x")) ~ / cL;|/(cd)| V s 

N r _ N r _ 

G(oj, x r , x)G(uj, x r , x') G(oj, x r , x)G[u>, x r , x") 

r—1 r—1 

x { g \G(", *, Vs)\ 2 E x '> Vs)G{^, x", y s ) - g \G{oj, x, y s )\ 2 G(u, x', y a )6(u>, x", y s )} 

s—1 s—1 s—1 

N r _ N r _ 

G(lo, x r , x)G(u), x r , x') G(uj, x r , x)G(u>, x r , x") 

r—1 r—1 

N s _ N s _ 

x l^G(u,x,y 3 )G(uj,x',y a ) ^ G(u, x, y s )G(uj, x", y s ) 

s=l s=l 
N s _ 

- ^2 G(uj, x, y s ) 2 G(u, x', y s )G(uj, x', y s )\ 



s=l 
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Here, T T = (J p 2 (t) di) and p T (t) is the probability density function of the random time delays. 
We will see that this quantity should be large so that the kernel T*T(x, x') is statistically stable. 
Using the method of Lagrange multipliers one can show for example that the maximal value of T T 
amongst all probability density functions p T compactly supported in [— T max ,T max ] is obtained for 
the uniform density over the interval and gives T T = 2T max . Since T max must be at most of the order 
of the recording time r max ~ T/2, to obtain a large value of T T one should take T large too (recall 
that in section 1.2 we had already assumed T to be large). 

Observe that in the above equation for the covariance, in each of the two terms on the right hand 
side we are summing N s terms with a minus sign and N 2 with a plus sign, and they are all of the 
same order. The contribution of the terms with a minus sign is therefore small, and we can use the 
results of section 1.3 to simplify the above equation into 

Cov(lC{x,x%fC{x,x")) ~ ^-L- J dwl/HlV 

x < sincf — \x — x'\ jsincf — \x — x"\ Jsincf — \x' — 
I \Cf) J Vc / Vc 

• 2 / W i /i\ . 2 / W i I, 

+ sine — \x — x \ ) sine — \x — x 
V c ) V c 

Similar computations for the case of stationary random sources, where the terms with a minus sign 
do not even appear, leads to 



x 



Cov(K,(x,x'),K,(x,x")) ~ — / dw|F(o;)| 2 w s 



<, sincf — \x — x'\ jsincf — \x — 
Vc / Vc 



"\)smc( -\x' - x"\ 
J Vc 



• 2 / w i n\ ■ 2 ( u i /, 

+ sine I — | x — x I sine I — | x — x 

V c / V c 

Let us continue the computations in the case of noise blended sources. Putting everything to- 
gether, we find that the standard deviation of the imaging functional at a point x is given by 



S(x) = 



1 



2ttT t 



with 



Ti(x, ui) — // sincf — \x — x'\] sincf — 
J Jb r L Vc o 1 Vc o 



x — x I I sine I — \x —x 

CO 



f — I sc - x'\)smc 2 ( — \x - x" 
V c / V c 



Sc- 2 (x')Sc- 2 (x")dx' Ax'' 



Again, we focus on the spatial integral Xi (and to simplify notations we drop the dependence on 

uj): the quantity we need to study is [^(x )/T r ] X ^ 2 for the noise blended sources and \T ■ X2(x)] 
for the stationary random sources. The computations presented below are performed in the noise 
blended sources setting. To obtain the corresponding standard deviation for stationary random 
sources it suffices to substitute the factor l/T T (or 1/a/TV) by T (or VT). For comparison with the 
average amplitude, recall that for stationary random sources the results obtained in the previous 
section have to be multiplied by a factor T. 
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We can rewrite the integral T 2 (x) as the sum of the two integrals 

J 1 (x)= [[ H 1 {x,x / ,x' r )5c- 2 {x')5c- 2 (x") dx' dx 

J J Br 

J 2 (x)= [[ H 2 {x,x' 1 x")5c- 2 {x')5c- 2 {x") dx' dx 

J J Br 

which is to say the double integral of the two kernels 
Hi(x,x',x") 
H 2 (x,x',x") 

applied to the perturbation (8c~ 2 {x')8c~ 2 (x")) , for every one of the three types of perturbation 
studied above. However, since 

J 2 (x)=ll(x), 

we will only need to analyze J\. 

For the following estimates it is important to keep separated the scale of the dimension of the 
perturbation (e) from the scale of the wavelength of the sources (77). Their relative amplitude will be 
specified, but to help the reader to keep track of the different orders, we stress that we will always 
have £ < n. 

3.1 Point singularities 

We return to the case of point perturbations introduced in the previous section and modeled by 
(2.1). In this case, even a very rude estimation is sufficient to obtain a bound which guarantees that 
this perturbation can be imaged. Since \sinc\ < 1, changing variables we get 

I J7i (je) I <£ 6 JJ |sinc(|aj — ex'|/?7)sinc(|a; — ex"\/rj) | dx' dx" 

<Ke 6 I sin^dx'le/^dsc', (3.2) 

JB 1 ( x /e) 

where the constant K = 4/3-7T comes from Cauchy-Schwarz inequality Since the sine function is 
bounded, we have obtained a bound of order e 6 . This is only a rough upper bound, but it is not 
necessary to look for an improvement since it is already of the same order of the integral of the 
second kernel, for which we have 

I 2 (0)=X 2 (0)~e 6 . 

Using the decay of the sine function, we can find near the perturbation {\x\ ~ e) a bound for T 2 of 
the same order. Oscillations are therefore of order e 3 / ' \/T\. 

Recall that the average value observed on the peak is of order e 3 for e <C n, so that the typical 
value observed remains of the same order due to the large factor T T . The same results holds true 
also for e ~ 77, namely the typical value observed is of the same order of the average value. 

Far from the perturbation the integrals of the two kernels decrease. Using the bound (3.2), the 
integral of the first kernel can be bounded for \x\ = 0{1) by 

\Ji(x)\ <Ke 6 [ n 2 \x - ex'\~ 2 dx' = K e 6 n 2 . 

JBx 



u 1 i\\ • ( I r>\\ ■ ( u 1 ' a 

— \x — x sine — \x — x smc — \x — x 

c / Vc / Vc 

• 2 S 10 1 t\\ ■ 2 ( w 1 
smc I — \x — x I smc I — \x — x 

VCn / Vcn 
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With some more work, one can show that this bound is sharp, at least for e <C n. This can be done 
using the Fourier representation of the sine function written in spherical coordinates u = u(r, 0, <f>) € 



sinc(A|a;|) 



1 
2 
1 

-in 



-iA£|as| 



-iAjajj cos(f?) 



sin(6>) dO 



o Jo 



sin(6»)d6id0 = 



4tt 



du 



(3.3) 



where S* 2 = {x £ M 3 : \x\ = 1} is the unitary sphere in R 3 . We can write 



Jx{x) = (4vr) 

=•6 



-3 



3 2 4tt 



s H(«-««')-V^"-« , ">*e < fl ( "'- , "" ),to d«dt;dto dx' dx h 



ii«-(u+«) du dvdw 



where we have used the assumption e <C rj. Simplifying this equation and using again (3.3), we get 
for \x\ = 0(1) 





2 e 6 r 


l i — x-u i 

/ e *> citt 




L Js 2 


~ 3 2 . 



47rsinc(|;c|/77) = 0(e 6 rf) 



As for the integral of the second kernel, the bound we have is of order e & rf. For e <C n we have 
therefore 



I 2 (x) = 0( £ V)- 

In the general case e < ?y, the above equation becomes an upper bound. 

Recall that the statistical average of the imaging operator far from the perturbation is of order 
e 3 n 2 . Assuming that T is large, but still 1 <C T T < l/rf, the above result implies that the typical 
value observed is at most of order e^n/y/T^ (it is exactly of this order for e <C rf). Therefore, the 
typical contrast is still at least of order y/T^/r}, allowing for a precise location of the perturbation 
(both when and e ~ rj). 



3.2 Line singularities 

Consider the case of line singularities, modeled by (2.2). Using rude estimations similar to the ones 
presented above, we could only bound the integral of the (absolute value of the) first kernel near the 
origin with something of order £ 4 ry 2 ln 2 (e). This means that, in order to be sure to able to image the 
perturbation, we would need to have e\ ln(e)| -C rj. But we can do better. 

Assuming simply e <C ??, we can approximate 

sinc^i-\/e 2 (.T' 2 + y' 2 ) + z' 2 ^j ~ sinc(|^'|/r?) = sinc(z'/?7) . (3.4) 
It is then possible to use the Fourier representation of the sine function 

sinc(Az) = 1. J e"^ d( = \ ^ e~^ z d( 
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to obtain the amplitude of oscillations. At x = we can rewrite the integral of the first kernel as an 
integral over C\ = B\ x [—1, 1] C R 2 x K, use (3.4) and integrate in x',y',x",y": 

Ji(0) = J J sinc(|a;'|/??)sinc(|a;"|/?7)sinc(|a; / -a;"|/77)d3;'da;" 

= £i ff sinc (^ V £2 ( x ' 2 + y' 2 ) + z ' 2 ) sinc (^ V £2 ( x " 2 + v" 2 ) + z " 2 ) x 



x si 

z+l 



(^^e 2 ((x> - x") 2 + {y' - y") 2 ) + (z' - z") 2 ) dx' dx' 



1 1 sine (z'/rj) sine (z"/r)) sine ((z' - z")/r]) dz'dz" , 



and using the Fourier representation introduced above we get 

Ji(0) = / f +l Iff e- 8C ' 2 '/" e - jC " 2 "/ r 'e- i;c ( z '- z ")/"dCdC'dC"dz'dz" 



^ sinc(i(C + C))sinc(-V - C)) dCdC'dC" 
e rj 1 I I sinc(Mi)dMi / sinc(ii2)du2 dC • 



n 2 r 1 /•(C+i)/'? /•(-C+i)/'? 

2 j-u( ( -i)/ n J(-C-i)/'? 
Since the function s i— > sinc(u) du is uniformly bounded in s, we get that 

Ji(0) = O(eV) ■ 

For the second kernel, we have 

^ 2 (o)=i 2 (o) = o( £ V). 

Remark that for |a;| ~ e we can still bound fluctuations in the same way, because (3.4) still holds, 
and the integral 

1/2 



(x - ex') 2 + (y- £y') 2 + (z - z') 2 ) dx' 



is maximal when x = 0. 

We see that fluctuations near a; = are of order Erj 2 /^/J\-. Since the average value observed 
at a; = for the imaging functional is of order Erf 2 , it is thanks to the large factor T T ^> 1 that 
we get the statistical stability of the operator. This means that the typical value observed on the 
perturbation remains of order erf 2 . 

When x is far from the perturbation, oscillations are even smaller. Indeed, we can proceed as in 
(3.2) to find a bound for the integral of the absolute value of the first kernel. We get 



Ji(x) < £ 



l 

4 



(- y/(x - £x') 2 + {y- £y') 2 + (z-z') 2 )|dx' dy' dz' 



sine 

2 

((x - £x'f + (y- £ V 'f + (z - zfY 1 ' 2 dx' 

-1 JB 1 



Since the integrand is bounded, the bound we get is of order e 4 r/ 2 . The second kernel is of a higher 
order, Ji(x) < 0(e i rj 4 ). This bound implies that the typical value observed far from the perturbation 
is of order at most E 2 r\l\jT\ 1 so that the contrast is at least of order r/y/T^/e. 
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3.3 Plane singularities 

We turn now to analyze fluctuations in the case of surface-type perturbations, modeled by (2.3). 
The difficult part is again to obtain good estimates on the integral of the first kernel, for which we 
use the Fourier representation of the sine function obtained in (3.3). At x = we have 



Ji(0) 



1 



(4tt) 3 
1 

(4tt) 3 

4e 2 
(4tt) 3 



;i [*'.(«+«)+*"•(*-•«>)] dM dv dw dx > dx i> 



a x '-( U+W ) dx , du 



(«+») dx > du 



dw 



dw , 



where the approximate equality holds for e <gc 77 and _L denotes the projection on the last two 
coordinates: x± = (y,z) <G M 2 . The integral in dx' ± is computed on (D £ )± = B\ € M 2 . We can 
rewrite also the integrals on S 2 as (twice the) integrals on the projection B\ £ M 2 , compute the 
integral in dx' ± and change variables: 



Ji(O) 



32e 2 
(4tt) 3 

2c 2 

7T 

2e 2 rf 



Bj 



Bi 



Bi 

Jl(|«_L 



yl — |i"±| 2 dw± 



w±\rj) 



Bi 



u>_l|/?7 



dtt 1 



1 dii; 1 



VT 



12 „2 



?7 2 ditj 



yl — \w±_\ 2 ri 2 dw± 



Here, Ji is the Bessel function of the first kind. Let us focus on the integral inside the square 
brackets. Observe that the origin of our system of coordinates is always inside B 1 / rj (w±). Changing 
to polar coordinates we have 



y 



■Mluxl) 



■\/l — \u± — w±\ 2 n 2 ditj 



2tt fP-u,(6) 

Ji(r)(f> w (0,r) drd6 



where we denote by <p w the square root term (written in polar coordinates) , and the function p w (9) 
takes values in [I/77 — \w±\, 1 / 77 + \w±\] C [0, 2/r]}. We claim that the integral term y is bounded. 
This can be proved integrating by parts in r. Remark that the square root term is concave (as a 
function of u±), take its maximum over the domain of integration B\ /^(itfjj at the center of the ball 
and is zero at the boundary. Therefore, for every fixed (6,w), the function <p w (0,-) is still concave 
and bounded by 1. Then, one easily obtains that the integral of the absolute value of its derivative 
in r is bounded by 2. Also, the antiderivative of Ji(r) is the Bessel function of order zero — Jo(r), 
which is bounded (the maximum of its absolute value is taken at r = 0, and Jo(0) = 1). Putting 
everything together, we get 



y 



2k 



< 



J o (0) + J o (0) 



r=0 
p w (9) 



J o (r)d r </) w (0,r)dr d9 



\d r <j) w {e,r)\dr dO < 6tt . 







This proves the claim. 
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Using again the boundedncss of the square root, we can bound the integral in dwi by ir/rj 2 . We 
have therefore obtained a bound for J7i(0) of order 0{e 2 rf). This is only an upper bound, but there 
is no need to look for an improvement, since it is already of a smaller order than the integral of the 
second kernel, for which we have 

J 2 (0)=X 2 (0) = O(e 2 r } A ]n 2 (e)). 

Therefore, 

I 2 (0)~eVln 2 (e) ■ 



Far from the perturbation, oscillations are even smaller. Denote x = (x,x±) € R x R 2 and 
u = (tii, u±) G R x R 2 ; the same notation is used for v and w. Let us look at the integral of the 
first kernel; following the computations presented above we have 

Ji(x) = J I [[[ e i7K«-»')-(«4i.)+(»-»")-{«-»)] dlidt , (lt0 da;' da;" 



S'2 



As 



2 



(4tt) 3 



32e 2 
(4tt) 3 



gi^[(as-a>') ±-(u+w) j_+(as-sc") ± -(v-w) J fa." 
S 2 J J Bi 

x ei L[ x {u +W ) 1+ x(v-w) ± ] dudvdw 

x e iHV 1 -l u -Ll 2 +V 1 -KI 2 ] 

x ^/l - |mj_| 2 ^/1 - |vj_| 2 ^l - |t/J^| 2 dttx duj. dwj_ 



so that 



F 2 6 /• 



Bi/,(i"j 



V ± yi — — w±\ 2 r) 2 — p — j-^dw 



Jb,,„(- wi ) «-L 



For | a; | ^> 1, the last two integrals above are now much smaller than the corresponding ones for 
a; = 0. This is due to the fact that for \x\ large, at least one of the two exponential terms, which 
have mean zero, is rapidly oscillating with respect to J\. We therefore have that J7i(a;) is at most 
of order e 2 rf. Far from the perturbation, the (sharp) bound we have on the integral of the second 
kernel is of the same order: Jy.{x) < 0{e 2 rf). We have obtained that 

T 2 (x)~e 2 r) 4 . 

Thanks to the large factor T T , fluctuations are therefore smaller than the average value given by the 
imaging functional, both on the perturbation and far from it. The typical contrast remains therefore 
of the same order as the average contrast, namely of order | ln(e)|. 
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4 Conclusions and comments 



We have analyzed the imaging functional given by (1.4)in the high frequency regime (rj <C 1) with 
respect to small perturbations (e 1) of the background velocity of propagation. Using a suitable 
disposition of the sources and receivers, we have been able to obtain quantitative estimates on 
the (average, with respect to the realization of the random time delays or the stationary random 
source signals) sensitivity of the imaging functional. The image presents a peak on the location 
of the perturbation, and the contrast is of order r\~ 2 for point perturbations, of order e~ l for line 
perturbations, and only of order | ln(e)| for surface perturbations. 

The most interesting result obtained in this paper concerns the quantitative analysis of the 
statistical stability of this functional, providing the typical contrast seen for the three perturbations 
considered. The question of stability of the imaging functional has been addressed in [DFGS12]: 
no quantitative analysis was carried out there, but it was shown that a condition for the statistical 
stability is that the quantities T (for stationary random sources) and T T (for noise blended sources) 
must be large. For random time delays uniformly distributed on the interval [— T maXl T maa: ], T T large 
means that T max must be large, which in turn implies that the recording time T ~ T max must be 
large. 

An important fact is that the typical contrast found only depends on the type of perturbation 
one is trying to image, and not on the method used. All results are described below for noise 
blended sources, but the corresponding contrast for stationary random sources are obtained simply 
substituting T T with T. 

We have shown that for point perturbations, both when e rj and e ~ rj, fluctuations due to 
the stochastic nature of the method are small, and the typical contrast is at least of order \fT r jjr] 
(y/T/rf for stationary random sources): point perturbations are easy to find. 

For line perturbation the situation is different. Wc can image with a satisfactory accuracy and 
reasonable recording time T> 1 only very thin line perturbations, e -C rj. The typical contrast in 
this case is at least of order \fTj-r\je. 

For plane perturbations the average contrast is quite poor, only of order | ln(e)|. However, for 
very thin perturbations , e <C rj, also the typical contrast is of the same order. 

These results are summarized in the following tables, where we list for the three type of perturba- 
tions considered the order of the average value given by the imaging functional and of the standard 
deviation at the center of the perturbation and far from it. 





(5c- 2 (x)) 
x = 


|sr| > 1 


S(x) 
x = 


\x\ > 1 


Points 
Lines 
Planes 


~ erj 2 
~ erj 2 |ln(e)| 


< e^rj 2 

< £ 2^2 

< erj 2 


ce 3 /vr; 

~ £TI 2 /VTt 

~er/ 2 \ln(e)\/VT; 


< eVVlV 

< s 2 ri/VTj 

< erj 2 



Table 1: Noise blended sources: mean (Sc~ 2 ) and standard deviation (S) of the estimated velocity 
perturbation at the center of the perturbation (x = 0) and far from it (|a;| ^> 1), in the regime 
e <C rj <C 1. The cases of point, line and disc singularities are displayed. 
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(6c- 2 (x)) 




5(x) 






x = 


\x\ > 1 


£C = 


|*| > 1 


Points 


~ Te 3 


< T£ 3 V 2 


~ VTs 3 


< \/Te 3 r/ 


Lines 


-Ten 2 


< ^2^2 


~ \/T£?7 2 


< VTe 2 77 


Planes 


~T£77 2 |ln(e)| 




~ VT£77 2 |ln( £ )| 





Table 2: Stationary random sources: mean (<5c~ 2 (a:)) and standard deviation (<S) of the estimated 
velocity perturbation at the center of the perturbation {x = 0) and far from it (|a;| 3> 1), in the 
regime s <C n <C 1. The cases of point, line and disc singularities are displayed. 
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